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ABSTRACT 

The Large Area Telescope (LAT) aboard the Fermi Gamma-ray Space Telescope provides an un- 
precedented opportunity to study gamma-ray blazars. To capitalize on this opportunity, beginning in 
late 2007, about a year before the start of LAT science operations, we began a large-scale, fast-cadence 
15 GHz radio monitoring program with the 40-m telescope at the Owens Valley Radio Observatory 
(OVRO). This program began with the 1158 northern (6 > —20°) sources from the Candidate Gamma- 
ray Blazar Survey (CGRaBS) and now encompasses over 1500 sources, each observed twice per week 
with about 4 mJy (minimum) and 3% (typical) uncertainty. Here, we describe this monitoring pro- 
gram and our methods, and present radio light curves from the first two years (2008 and 2009). As 
a first application, we combine these data with a novel measure of light curve variability amplitude, 
the intrinsic modulation index, through a likelihood analysis to examine the variability properties of 
subpopulations of our sample. We demonstrate that, with high significance (7-cr), gamma-ray- loud 
blazars detected by the LAT during its first 11 months of operation vary with about a factor of two 
greater amplitude than do the gamma-ray quiet blazars in our sample. We also find a significant 
(3-cr) difference between variability amplitude in BL Lacertae objects and flat-spectrum radio quasars 
(FSRQs), with the former exhibiting larger variability amplitudes. Finally, low-redshift (z < 1) FS- 
RQs are found to vary more strongly than high-redshift FSRQs, with 3-cr significance. These findings 
represent an important step toward understanding why some blazars emit gamma-rays while others, 
with apparently similar properties, remain silent. 

Subject headings: BL Lacertae objects: general — galaxies: active — methods: statistical — quasars: 
general — radio continuum: galaxies 
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1. INTRODUCTION 

The rotating super-massive black holes that power ac- 
tive galactic nuclei (AGN) somehow accomplish the re- 
markable feat of channeling energy derived from their 
rotation and accretion disks into two relativistic jets 
oppositely-directed along the spin axis. In spite of in- 
tensive observational efforts over the last four decades, 
the detailed mechanism of this process has remained 
elusive, and, although several processes have been sug- 
gested, we are still largely ignorant of the composition 
of the jets and the forces that collimate them. The 
first detailed collimation mech anism to be proposed wa s 
that of a "de Laval" nozzle (jBlandford fe Reesl Il974h . 
which is now known to be a likely cause of re-collimation 
on kpc scales, but not of the initial collimation, 
which, as revealed by Very Long Baseline Interferomc- 
try (VLBI), clearly occurs on sub-parsec scales. Other 
early theories, which involve magneto-hydrodynamic 
winds (Blandf ord fe Znaieklll977l ) an d /or magnetic fields 
threa ding the inner accretion disk (jBlandfor d fe Payne! 
1982), remain the most promising approaches to a full 
understanding of the phenomenon. 

An observational difficulty is that, except in a few 
cases (e.g. M87), radio observations, which provide 
the most detailed images of active galaxies, only probe 
the relativistic jets down to the point at which the 
jets become optically thick at a point some light-weeks 
or light-months from the site of the original collima- 
tion. Higher-frequency observations are needed to probe 
deeper into the jets, although interstellar scintillation 
observations do in some cases reveal the presence of 
radio emission features in so me AGN that are ~ 5- 
50 m i cro-arcseconds in exten t (Ked ziora-Chudczer et al 
1997t iDennett-Thqrpe fc de Bruvnl 120001: iJauncev et aL 
200q |Rickett et alj|2002l 120061: lLovell et al.ll2008l) . which 
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can be very persistent (jMacauart fc de Bruvnl 120071 ). 
These mysterious, very high brightness temperature fea- 
tures are by no means understood, and are certainly 
of great interest. At optical wavelengths, rapid swings 
in the polarization position angle have been used to 
tie together flux density variations at TeV energies 
and y ariations at millimeter wavelengths (Marschc r et al.l 
[2008h . At very high energies of hundreds of GeV 
to TeV, very rapid variations down to timescales of 
minutes have been observed b y the HESS, MAGIC 
and VERITAS instruments (e.g. lAharonian et al.1 [20071 
Aharonian et al.l [20091: lAcciari et al.l 120091: lAcciari et al.l 



20101) . Full three-dimensional (non-axisymmetric) mag- 



netohydrodynamic relativistic simulations are now being 
carried out that enable detailed interpretation of the ob- 
servations over the whole electromagnetic spectr um (e.g. 
IMcKinnev fc Blandford! 120091 IPenna et al.ll2010D . 

Relativistic beaming introduces complications in ob- 
servational studies of relativistic jets. The continuum 
emission is strongly beamed along the jet axis, introduc- 
ing strong observational selection effects. Those objects 
having jets that are aligned at a small angle to the line 
of sight are collectively known as "blazars." Small vari- 
ations in the angle between the jet axis and the line of 
sight result in a large range of observed properties, such 
as apparent luminosity, variability, and energy spectrum. 
Strong boosting of the continuum synchrotron emission 
from the jet also frequently swamps optical line emission, 



making it difficult or even impossible to obtain a redshift 
for the source. As a result, blazars are subdivided into 
two classes: flat-spectrum radio quasars (FSRQs) and 
BL Lacertae objects (BL Lacs). The former class con- 
tains blazars dominated by strong broad emission lines 
while the latter class contains those blazars with spectra 
dominated by their continuum emission, and hence weak, 
if any, emission lines and very weak absorption lines, or 
no lines at all. The large variations in the energy spec- 
trum make it difficult to study many blazars over the 
whole electromagnetic spectrum. As a result the study 
of large, carefully-selected samples is necessary to deter- 
mine the physical processes and conditions of the parent 
population. As relativistically boosted emission can be 
detected even from high-redshift sources, any intrinsic 
scatter in jet properties and scatter due to relativistic 
beaming is additionally convolved with any cosmological 
evolution of the black holes giving rise to the jets and 
their environment. It is therefore not surprising that the 
study of the population properties of relativistic jets has, 
to this day, been sparse at best. 

The launch of the Fermi Gamma-ray Space Telescope 
in June of 2008 provides an unprecedente d opportunity 
for t he systematic study of blazar jets (|Atwood et alJ 
2009). Its Large Area Telescope (LAT) observes the 
sky at energies between 100 MeV and a few hun- 
dred GeV. In this energy range relativistic particles 
can be probed through their inverse Com pton emission 
in the case of ele ctron/positron jets (e.g. iDermer et al.1 
H9<&ISlkora et alJll994b Blandford fc LevinsorJll995[ ). or 
a combination of pionic emission from primaries and in- 
verse Compton emission from cascade-produ ced leptonic 
secon daries in the case of hadronic jets (e.g. iMannheiml 

mil). 

Blazars comprise the most numerous class of ex- 
tragalactic GeV sources associated with lower-energy 
counterparts: the first- year Fermi point source cata- 
log (1FGL) contains 1451 sources, of which 596 have 
been ass ociated with blazar s in the first AGN catalog 
(1LAC) (|Abdo et al.l I2010al lbh. As the LAT completes 
one survey of the whole sky every 3 hours, it can provide 
continuous monitoring of all gamma-ray bright blazars, 
although with variable cadence that depends on the inte- 
gration time necessary to detect each object (which can 
range from a single satellite pass to many months, de- 
pending on the average flux density of the object and its 
activity state). 

The exact location of the gamma-ray emission re- 
gion and its proximity to the central black hole re- 
main subjects of debate. Two possible models of 
the GeV emission region are that this emission comes 
from a "gamma-sphere" clo se to the base of the jet 
(Blandfor d fc LevinsonllT"995[ ). or that it comes from the 
same shocked regions that are responsible for the radio 
emissio n seen in VLB! obse rvations much further out in 
the jet (|Jorstad et a l. 2001). If the former model is cor- 
rect then the gamma-ray observations might well provide 
evidence of the initial collimation mechanism. 

The testing of models of the location, structure, and 
radiative properties of the gamma-ray emission region in 
blazars requires, in addition to the Fermi observations, 
supporting broadband observations of likely gamma-ray 
sources in various activity states. Such multiwavelength 
efforts can occur in two modes: 
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1. regular monitoring of a preselected, statistically 
complete sample of likely gamma-ray-bright ob- 
jects, independent of their gamma-ray activity 
state; and 

2. intensive observations of archetypal objects or ob- 
jects exhibiting unusual behavior 

The blazar monitoring program we discuss here is fo- 
cused on the first mode. In anticipation of the unique 
opportunities offered by the Fermi-hAT sky monitor- 
ing at gamma-ray energies, three years ago we be- 
gan the bi-weekly monitoring of a large sample in- 
cluding 1158 likely gamma-ray loud blazars, prese- 
lected according to uniform criteria, with the Owens 
Valley Radio Observatory (OVRO) 40-m telescope at 
15 GHz. We also apply our observations in studies of 
the second mode through Fermi-hA T multiwave l ength 
campaigns for flaring sources (e.g. JAbdo et al.l 120091 : 
\Fermi-hAT Collaboration et al]|2010f l) and through col- 
laboration with the F-GAMMA project, a complemen- 
tary effort representing the second mode, focused on ra- 
dio and sub- m m spectral monit o ring o f about 60 promi- 
nent sources (jAn gclakis e t al.l 120101 : iFuhrmann et al.l 
12001 . 

The sample that we are studying with the OVRO 40- 
m telescope is statistically well-defined and large enough 
to allow for statistical analyses and comparisons of sub- 
samples. In addition, as the 40-m telescope is dedicated 
full time to this project, the cadence is high enough to 
allow sampling of the radio light curves on timescales 
comparable with those typically achieved by Fermi-hAT 
for bright gamma-ray blazars and in this sense the 40-m 
and Fermi-hAT are ideally matched. The combination 
of sample size and cadence is unprecedented, making this 
by far the largest monitoring survey of radio sources that 
has been undertaken to the date of writing. 

Data from this program, in combination with Fermi 
observations, will allow us systematically to derive the 
radio and radio/gamma ray observational properties of 
the blazar population, including 

• the radio variability properties of the blazar popu- 
lation, their dependence on redshift, spectral clas- 
sification, luminosity, and gamma-ray activity; 

• any differences between the radio properties of 
gamma-ray loud blazars and blazars with similar 
radio luminosity which have not been detected by 
Fermi ; 

• the properties (e.g. significance of correlation and 
the length and sign of any time delays) of cross- 
correlations between radio and gamma-ray flares 
of gamma-ray loud blazars; and 

• the combination of radio properties, if one exists, 
that can predict the apparent gamma-ray luminos- 
ity of a blazar (which, in turn, could be used to 
derive blazar gamma-ray luminosity functions from 
radio luminosity functions). 

Such a systematic study of radio and radio/gamma pop- 
ulation properties should allow us to address a series 
of long-standing questions on the physical properties of 



blazar jets, including the location, structure, and radia- 
tive properties of the gamma-ray emission region, and 
the collimation, composition, particle acceleration, and 
emission mechanisms in blazar jets. 

In this paper we describe in detail the 40-m telescope 
15 GHz monitoring program, we present results from the 
first two years of the program (2008 and 2009), and we 
derive the variability amplitude properties of the blazar 
population at 15 GHz. Studies of other blazar population 
radio and radio/gamma pro perties will be discussed in 
upcoming publications (e.g. iMax-Moerbeck et al.l 120111: 
Pavlidou etaLll2011l lAbdo et all 1201 H IFuhrmann et al.l 



2011). 



The remainder of this paper is organized as follows. In 
§ 2, we discuss the telescope and receiver and our mea- 
surement procedures. In § 3 we discuss the method of 
operation. In § 4 we discuss our sample of sources and 
observing strategy. In § 5 we discuss data editing and cal- 
ibration. Our results, including light curves, the deriva- 
tion of variability amplitudes for the blazar population, 
and population studies using this analysis are presented 
in § 6. We summarize and discuss our conclusions in § 7. 

2. TELESCOPE AND RECEIVER 
2.1. Optics 

The OVRO "40-m" telescope is actually a 130-foot- 
diameter f/0.4 parabolic reflector with approximately 
1.1 mm rms surface accuracy on an altitude-azimuth 
mount. A cooled receiver with two symmetric off- axis 
corrugated horn feeds is installed at the prime focus. The 
telescope and receiver combination produces a pair of ap- 
proximately Gaussian beams (157" FWHM), separated 
in azimuth by 12.'95. We refer to these two beams, some- 
what arbitrarily, as the "antenna" beam and the "refer- 
ence" beam, or ant and ref. The receiver selects left-hand 
circular polarization, so linearly polarized sources of all 
orientations may be monitored in total intensity. By ob- 
serving compact sources of known flux density, we find 
the aperture efficiency, t]a ~ 0.25. 

This relatively low aperture efficiency is due to deliber- 
ate underillumination of the dish by the feed — for mon- 
itoring observations of a large sample of objects aiming 
at flux density measurements repeatable to within a few 
percent we must consider the trade-off between aperture 
efficiency and pointing accuracy. Underillumination of 
the antenna increases the beamwidth and reduces sus- 
ceptibility to pointing errors relative to more fully illu- 
minating the antenna, in addition to reducing exposure 
to thermal noise from ground spillover. Experience has 
shown that we are operating at close to the optimum 
illumination for the most efficient use of the telescope: 
increasing the aperture efficiency gains little because the 
thermal noise is already acceptably low for observing the 
objects in our monitoring sample. The on-source duty 
cycle is about 20% — a factor of two in efficiency is lost 
due to Dicke switching, and the rest of the lost time is 
due to slewing and calibration. 

When the 40-m telescope moves in elevation, gravity 
deforms its surface, changing the antenna gain and focus 
location. The entire feed/receiver system can be moved 
along the optical axis to adjust the focus. The opti- 
mum focus position as a function of elevation is measured 
about once per year, but has not been found to vary sig- 
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Fig. 1. — Block diagram of the Ku-band receiver. 



nificantly except when the receiver has been removed and 
reinstalled during maintenance. Due to thermal effects, 
the optimum focus also varies slightly between day and 
night operation and with the angle between the telescope 
structure and the Sun. 

2.2. Receiver 

A block diagram of the receiver is shown in Figure [TJ 
The receiver operates in the Ku band with a center fre- 
quency of 15.0 GHz, a 3.0 GHz bandwidth, and a noise- 
equivalent reception bandwidth of 2.5 GHz. The receiver 
noise temperature is about 30 K, and the typical sys- 
tem noise temperature including CMB, atmospheric, and 
ground contributions is about 55 K. 

In order to make the most efficient use of the telescope, 
a Dicke-switched dual-beam system is used. A ferrite 
switch alternately selects between the ant and ref beams 
and delivers the difference between the two, which is the 
switched power, i.e. the difference between the power in 
ant beam and the power in the ref beam. We designate 
this power difference by £, i.e. 



Pre 



(1) 



Although Dicke-switching halves the time spent observ- 
ing the object it is more efficient than using a single-beam 
and scanning the telescope across the source because the 
integrated signal from the source is higher and hence flux 
densities can be measured faster, and in addition, as de- 
scribed below, Dicke-switching removes large systematic 
errors. 

The receiver front end consists of a cooled (T^80 K), 
low-loss ferrite RF Dicke switch followed by a cryogenic 
(T~13 K) HEMT low-noise amplifier. This is followed by 
additional room-temperature amplifiers, a 13.5-16.5 GHz 
band definition filter, and an electronically-controlled at- 
tenuator used to adjust the overall gain of the receiver. 
The signal is detected directly using a square law detec- 
tor diode. The detected signal is digitized with a 16-bit 
analog-to-digital converter and then recorded. 

From 2007 until 2008 November, several receiver cal- 
ibrations were performed. Beginning in 2008 Novem- 
ber, approximately monthly calibrations were performed 
to monitor receiver performance. These calibrations in- 
cluded Y-factor measurements to characterize receiver 
temperature, sky dips to measure atmospheric optical 
depth and to determine the ground spillover, calibration 
diode effective temperature measurements, and observa- 
tions of calibration sources to measure the aperture effi- 
ciency. 

2.2.1. Measurement procedures 



In typical radiometry observations on the 40-m tele- 
scope we use three procedures: (i) flux density measure- 
ments, (ii) measurements of a calibration noise source, 
and (iii) pointing measurements on a nearby bright 
source. The receiver output voltage is integrated and dig- 
itally sampled at 1 kHz, synchronously with the 500 Hz 
Dicke switching rate. Alternate millisecond samples are 
subtracted to demodulate the Dicke switching, and the 
results are accumulated into 1-second averages. In addi- 
tion to the demodulated outputs, the sum of the alter- 
nate samples, i.e. the total powers in both the ant and 
ref beams are also recorded. 

2.2.2. Calibration diode 

A pair of calibrated noise diodes, referred to as the 
"noise" and "calibration" diodes, are connected to the 
main beam input via directional couplers to the Dicke 
switch. These noise diodes provide an excess noise ratio 
of (31 ± 1) dB from 12-18 GHz with stability of about 
0.001 dB/K. The diodes provide two calibration levels — 
one of power comparable to the system temperature and 
one attenuated to provide power comparable to the astro- 
nomical sources we are observing. The equivalent noise 
temperatures of the noise and calibration diodes at the 
receiver input are about 67 K and 1 K, respectively, sta- 
ble to < 1%. 

2.3. Pointing 

The 40-m telescope is equipped with encoders on the 
azimuth and elevation shafts and with two orthogonal 
tilt meters located in the teepee of the telescope in the 
alidade above the azimuth bearing. These four sets of 
readings are combined in a pointing model that generates 
encoder azimuth and zenith angle offsets based on the 
requested position on the sky. The pointing model has 
9 terms for the azimuth angle correction and 5 terms for 
the zenith angle correction, 
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'model = M sin 9 + A 2 + A 3 sin <f> cos 9 

+ A 4 cos <j) cos 9 + A 5 cos + A 6 sin <f> sin 9 (2) 
+ A 7 cos 4> sin 6 + A$ sin (40) + A g T LR cos 9 

A9 mode i = Zi + Z 2 sin 9-Z 3 cos (fi+Zi sin 4>+Z 5 T A f (3) 

Here, <fi an d 9 are the requested azimuth and zenith an- 
gles, A(j) mode i and A9 model are the pointing model cor- 
rections for the azimuth and zenith angles, Ai and Zi are 
the pointing model coefficients, and Taf an d Tlr are the 
aft-forward and left-right tilt meter readings. 

We have found that the pointing model terms drift 
slowly with time. Figure [2] shows the residual offset be- 
tween the pointing model and the actual requested posi- 
tion for 2008 and 2009. The sharp steps in the average 
offset correspond to adjustments in the pointing model. 
We adjust the pointing model two to three times per 
year to minimize the scatter in the offset and maintain 
an average offset less than about 0.'5 to ensure accurate 
pointing. Early in 2008 and at the end of 2009, the offset 
approached 1', but because the scatter did not increase, 
there was no substantial impact on data quality. 

In addition to the pointing model correction, at least 
once per hour we measure the pointing offset between 
a bright pointing calibrator and the model prediction. 
This measures the effect of wind and thermal loading. 
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100 200 300 400 500 600 
Days since MJD 54466 (Jan 1, 2008) 

Fig. 2. — Residual error between the pointing model and the 
actual requested position, plotted in week-long bins. The plotted 
data and errors are the weekly means and standard deviations of 
the pointing offsets measured by the pointing calibrations. 

In early 2009, we determined that these pointing offsets 
have the accuracy we require only at separations up to 
about 30° from the position where the pointing offset 
was measured. Because of this effect, after MJD 54906 
(2009 March 16), care was taken when scheduling to en- 
sure that flux density measurements were always made at 
separations of less than 15° from the pointing offset mea- 
surement. Prior to this, no such limit was in place. We 
have discarded flux densities measured with a separation 
of more than 30°. 

The pointing offset is measured by performing 3-point 
cross-scans of the calibrator in both azimuth and zenith 
angle and fitting a fixed-width Gaussian beam profile to 
each axis to determine the position of the peak. A point- 
ing offset measurement is considered invalid if its signal- 
to- noise ratio is less than 2, or if the offset indicates 
that the peak was outside the span of the cross-scan, 
± FWHM/2. Several iterations are attempted, moving 
the cross-scan center by up to FWHM/2 after each at- 
tempt, allowing offsets less than the FWHM (157") to 
be measured reliably. 

3. METHOD OF OBSERVING 
3.1. Double Switching 

The observing method we used follows closely 
tha t which we de v eloped and described in detail 
in iReadhead et all (fl~989 ) . and also discussed in 
lAngelakis et al.l 1 20091V To remove the large varying 
total power signal and minimize the effects of the at- 
mospheric fluctuations, ground spillover, and gain fluc- 
tuations we use a "double switching" approach, with a 
Dicke switch operating at 500 Hz, and azimuth switch- 
ing in which we alternate the beams on the object of 
interest. The advantage of double switching is that large 
variable signals are eliminated and the disadvantage is 
the loss of a factor of two in sensitivity — a factor of y2 
lost through observing the object only half the time, and 
another factor of v2 lost through the noise introduced 
by subtracting off the reference field. 

3.1.1. Dicke Switching 

The most important benefit of Dicke switching is the 
removal of the large, slowly varying total power signal, 
which is made up of contributions from ground, atmo- 



sphere, and receiver thermal noise. Variations in the gain 
of the low noise amplifier cause variations in the large to- 
tal power signal, and in addition the signals themselves 
vary slowly with time and with the position of the tele- 
scope. The resulting large variations in power limit the 
sensitivity of the receiving system. Ground spillover, like 
gain variations, contributes directly to the system noise, 
but the effect is difficult to quantify due to the com- 
plexity of the far sidelobes of the telescope beam. Dicke 
switching removes these large slowly-varying signals. 

A second benefit of Dicke switching is the reduction 
of noise due to the rapidly varying atmosphere above 
the telescope. With a beam separation of 12.'95, and 
for a water vapor scale height of 1.5 km, 75% of the 
total mass of water vapor seen by the telescope lies in 
the overlapping portions of the two beams. This fraction 
does not change substantially with scale height, dropping 
only to 72% (70%) for a water vapor scale height of 2 km 
(2.5 km). So Dicke switching reduces the effects of the 
varying atmosphere by about a factor of 4. 

A third benefit of Dicke switching is the relatively short 
observing time compared to the time required to scan 
across a radio source with a single be am. A detailed 
discus sion of these benefits is given in IReadhea d et al.l 
([198^ . 

3.1.2. Beam Switching and Flux Density Measurements 

While Dicke switching does much to reduce the large 
error terms due to the atmosphere, the ground, and gain 
fluctuations in the receiver, it does not remove linear 
drifts in any of these quantities and the situation can be 
further improved by beam switching. Beam switching in 
azimuth is optimum because by maintaining a constant 
elevation we minimize changes to the atmospheric and 
ground spillover signals and thereby maximize their can- 
cellation. We therefore adopt the same "doub le switch- 
ing" technique used by IReadhead et ail (H989), in which 
we alternate the two beams on the object of interest, 
and hence remove both the constant term and any linear 
drifts in the power from these unwanted components of 
the signal. 

The procedure we use for measuring fl ux densities is 
identi cal to that described in detail in IReadhead et al.l 
(1989) so we do not repeat all the details here, but give 
only a summary. To begin with the ref beam is posi- 
tioned on the source for 8 seconds, and the power differ- 
ence, !;a, is recorded. Then the ant beam is positioned 
on the source for 8 seconds and the power difference, £b , 
is recorded. With the ant beam still on the source a 
second 8 second observation is then made and the dif- 
ference, is recorded. Finally the ref beam is again 
positioned on the source for a final 8 second period and 
the difference, £d is recorded. Thus we spend a total 
of 32 seconds actually integrating on the source for each 
flux density measurement. Of course, slewing and set- 
tling times have to be allowed for at the beginning of the 
A, B and D integrations, so that the total time required 
for the flux density measurement is about one minute. 

The corresponding flux density is given by 



Sib = j{Cb + £c - £a - S,d) 



(4) 



where k is the calibration factor required to turn digitizer 
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units into Jy, and the rms error is given by 

(715 = Jv^4 + CT B + CT C + Cr D 

The calibration factor consists of a relative calibration 
factor that is computed for each measurement (§ 15.21) 
and an absolute calibration factor (§ I5.3[) . 

The four measurements also contain interesting infor- 
mation on the stability of the instrument and, more im- 
portantly, the atmosphere, during the observations. For 
each flux density measurement, we therefore also com- 
pute two other quantities — one that we call the "switched 
power," if>, given by 

^ = 1^B+^D-U-^C) (6) 

and the other that we call the "switched difference," //, 
given by 

ii= \{tc + tn-U-tB). (7) 

Both ip and /i should be zero in the absence of gain or 
atmospheric drifts so we use these as a way of estimating 
such variations in our error model (§ 15 .4[) and to reject 
badly-contaminated measurements f§ 15.1. 5[) . The uncer- 
tainties in if) and fj, are clearly given by equation (|S|). 

3.2. Confusion 

For sources at galactic latitude \b\ > 10° most of the 
reference fields are empty, but there are some objects that 
are contaminated by confusion introduced by other radio 
sources in the field. Fortunately since we are observing 
bright sources confusi on is not a problem. At 15.2 GHz, 
iWaldram et al.l (|2010f ) report a differential source count 
n(S) ~ 51(Syjy)~ 2 - 15 Jy -1 sr _1 with no deviation to a 
completeness limit of 5.5 mJy. Assuming that the effect 
of source clustering is negligible, the expected number 
of confusing sources detected at or above a flux density 
limit S c in either the ant or ref beam is 

N(S C ) = / n(S)n(S)dS (8) 

where Q(S) is the beam solid angle with antenna gain 
sufficient to detect a source of flux density S at the S c 
level. For a beam-switched flux density measurement, 
the expected number of confusing sources in the main 
or either reference beam is then v = N(S C ) + 2N (S c i/2) 
where the confusion limit is higher in the reference beams 
results because each is integrated only half as long as 
the main beam. Considering the confusing sources to be 
independently distributed among the observed fields via 
a Poisson process, the probability that a beam-switched 
flux density measurement includes one or more confusing 
sources is 1 — e~ v . 

Table [T] shows the probability of a confusing radio 
source lying in the main field or either reference field 
of a single flux density measurement, as well as the ex- 
pected number of contaminated sources in our 1158 ob- 
ject sample. Here, we have treated the ant and ref beams 
as identical 157" FWHM Gaussian beams and neglected 
reference field rotations with parallactic angle. The lat- 
ter approximation is justified because we observe sources 
at approximately the same local sidereal time each cy- 
cle, limiting the parallactic angles at which sources are 



TABLE 1 

Confusion 



Flux Density Limit (mJy) 


Probability 


Sources Affected a 


100 


8.4 X 10" 4 


1 


50 


1.9 x 10~ 3 


2 


20 


5.3 x 10~ 3 


6 


10 


1.2 x 10~ 2 


14 



a Expected number of contaminated program sources, consider- 
ing a source contaminated if a confusing source is found in the 
source field or either of its two reference fields. 



observed. Because only about 1.2% of our sources are 
likely to be contaminated even at 10 mJy level (~ 3% 
of the median flux density of sources in our sample) , we 
may safely ignore the effects of confusion in our statisti- 
cal analyses. 

4. OBSERVATIONS 
4.1. Source Selection 

The selection of the core sample for our monitoring 
program was driven by three considerations. First, since 
we are interested in the detailed study of the radio vari- 
ability properties of the blazar population and the de- 
pendence of these properties on other observables such 
as redshift, the sample should be sufficiently large to al- 
low division in subsamples (e.g. in redshift or luminosity 
bins) with enough members to derive confidently the sta- 
tistical properties in each. 

Second, to allow for the evaluation of the confidence 
level of any correlations or variable dependencies identi- 
fied in our data through Monte-Carlo simulations, and 
the generalization of our findings to the blazar popula- 
tion, the sample should be well-defined statistically, using 
uniform and easily reproducible criteria. 

Finally, one of the major goals of our monitoring pro- 
gram is the cross-correlation of 15 GHz light curves with 
Fermi gamma-ray light curves. For this reason we would 
like our sample to include a large number of gamma- 
ray loud-blazars. On the other hand, we would also like 
to be able to address the question of why some blazars 
are gamma-ray loud while other blazars, with apparently 
similar properties, are not. For this reason we would like 
our sample to be preselected — before Fermi data bias 
our understanding of what constitutes a likely gamma- 
ray-loud blazar — and, ideally, to include a comparable 
number of blazars which are not gamma-ray loud. 

Blazars in the Candidate Gamma-Ray Blazar 
Survey (CGRaBS) satisfy all of the requirements 
above (jHealev et al.l l2008h . CGRaBS blazars were 
selected from a flat-spectrum parent sample (complete 
to 65 mJy flux density at 4.8 GHz and radio spectral 
index a > —0.5 where S cx v a ) by a well-defined figure- 
of-merit criterion based on radio spectral index, 8.4 GHz 
radio flux density, and X-ray flux based on counts/s in 
the ROSAT All Sky Survey, to resemble blazars that 
were detected by the Energetic Gamma-Ray Experiment 
Telescope (EGRET, the precursor of Fermi-L AT). The 
CGRaBS sample is a total of 1625 active galactic nuclei 
(AGN) over the whole sky outside a ±10° band around 
the galactic plane. This sample was compiled before 
the launch of the Fermi and was expected to contain a 
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large fraction of the extragalactic sources that would be 
detected by Fermi-L AT . 

The core sample for our monitoring program consists 
of the 1158 CGRaBS sources north of declination —20°. 
As published, our subset of the CGRaBS sample con- 
tains 812 FSRQs, 111 BL Lacs, 53 radio galaxies, and 
182 objects without spectroscopic identification. In our 
analysis we use redshifts from the CGRaBS publication, 
which covered 81% of the sample (100% of FSRQs, 49% 
of BL Lacs). Recent spectroscopy has improved the com- 
pleteness of the sample to 886 FSRQs, 122 BL Lacs, 60 
radio galaxies, and 88 objects without spectroscopic iden- 
tification, with redshifts now available for 87.5% of the 
sample (100% of FSRQs; 53% of BL Lacs). The median 
15 GHz flux density for sources in our sample ranged 
from about 20 mJy to 30 Jy with a median of 325 mJy 
during the observation period described in this paper. 

In 1L AC, 709 AGN w ere associated with 1FGL 
sources (jAbdo et al.l l2010bh . Although continued im- 
provements in evaluating the probability of radio coun- 
terpart identification have caused some source associ- 
ations to vary in estimated significance and continued 
optical observations have improved the completeness of 
the typing and redshifts, we adopt here the id entifica- 
tions published in the 1LAC (|Abdo et al.ll2010bft . These 
identifications include 291 CGRaBS sources (221 of our 
subset) that were assoc iated with sources detected in 
gamma-rays by Fermi (Abdo et al. 2010b). Of these, 
263 (199 of our subset) were considered "clean" associa- 
tions, meaning only one source was associated with the 
gamma-ray source and the association probability was 
greater than 80%. CGRaBS sources made up 44% of 
the clean associations in the first-year Fermi AGN cata- 
log. This number is thus far smaller than anticipated; in 
the 11-month 1LAC sample only ~ 16% of the CGRaBS 
sources were detected and a large number of blazars not 
in CGRaBS have been detected. This suggests that the 
CGRaBS (EGRET-like) blazar sample is substantially 
different from that seen in the early Fermi mission. This 
finding represents a unique opportunity to investigate 
why gamma-ray activity is found only in certain blazars, 
and for this reason we retain in our monitoring program 
all of the blazars in our original core sample even if they 
have not yet been detected by the LAT. However, in 
order to optimize the potential for studies of the cross- 
correlation between radio and gamma-ray light curves, 
we have since added (and we continue to add) to our 
monitoring program all new Fermi-LAT blazars north of 
—20° declination that are not CGRaBS members. 

Several bright, stable non-blazar sources are included 
in our program to provide flux density calibration and 
to monitor instrumental variability. These are 3C 48, 
3C 161, 3C 286, DR 21, and NGC 7027. In addition to 
the stable sources, a number of bright sources are used to 
calibrate pointing. These sources need not exhibit stable 
flux density, but need be brighter than about 100 mJy 
to permit pointing offsets to be measured reliably. 

In addition to the core samples of blazars discussed 
above we have added further small samples of objects 
to our bi-weekly monitoring program, including (i) any 
objects not already included in our sample that are being 
studied in the F-GAMMA or VERITAS programs; (ii) 
a variety of galactic objects, such as microquasars and 
cataclysmic variables; and (iii) a few bright radio galaxies 



that show interesting jet properties. We are continually 
adding sources of interest to our monitoring sample, so 
that by now the sample comprises over 1500 objects that 
are monitored twice- weekly. 

4.2. Scheduling 

The large number of sources being observed requires 
the development of strategies to optimize the use of the 
telescope and minimize the effect of known systematic er- 
rors. The principal systematic errors we try to minimize 
are gain variations, atmospheric optical depth variations, 
and pointing errors. To achieve this optimization while 
minimizing slew times and dead times between observa- 
tions requires careful planning. Due to the size of our 
sample the scheduling must be automated. 

Schedules are arranged to ensure that sources are ob- 
served between zenith angles of 20° and 60° whenever 
possible. This is done for a number of reasons: 

1. the figure of the telescope was set for maximum 
gain in this elevation range; 

2. at zenith angles less than about 20° the telescope 
has to move rapidly to track an object and pointing 
accuracy can be degraded; 

3. at zenith angles greater than about 60° ground 
spillover increases significantly with decreasing el- 
evation; 

4. it is desirable to minimize the variation in atmo- 
spheric optical depth on our sources so as to mini- 
mize this particular source of error; and 

5. we try to minimize telescope slew times by observ- 
ing to the south and east in a limited elevation 
range. 

In the scheme we have developed, the sky is divided 
into 192 cells, each with a diamet er < 20°, using th e 
HEALPix mesh with N side = 4 (jGorski et al.1 12005[ ). 
Each source is assigned to a cell. From the sources in 
each cell, a pointing calibrator is selected using the fol- 
lowing criteria, applied in order: (1) if there is a flux 
calibrator in the region, this source is selected; (2) if one 
or more sources in the region have a flux density larger 
than 500 mJy, the one which minimizes the average an- 
gular distance to all the sources in that region is selected; 
(3) the source with the largest flux density in the region 
is selected. For these flux density comparisons, the me- 
dian flux density of the source during the previous year's 
observations is used. 

Sources within the region are ordered to minimize slew 
time, using a direct search to find the optimal order for 
regions with fewer than 9 sources and simulated anneal- 
ing for regions with 9 or more sources. A second opti- 
mization step determines the order in which the regions 
are scheduled using a heuristic algorithm in which re- 
gions are observed within a fixed zenith angle range and 
regions to the south have priority. The total sample is 
observed in three days. 

Prior to MJD 54906 (2009 March 16), a scheduling al- 
gorithm was used that did not enforce an angular sepa- 
ration limit between pointing calibrators and subsequent 
flux density measurements. 
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5. DATA EDITING AND CALIBRATION 

5.1. Data Editing and Flagging 

Editing and removal of corrupted data is performed 
using both automated and manual filters. 

5.1.1. Wind, Sun, Moon, and Zenith Angle Cuts 

Under high winds there is a systematic reduction 
in observed flux densities due to mis-pointing and 
poor tracking. Observations when the wind speed ex- 
ceeds 6.7 m • s^ 1 (15 mph) are discarded. To pro- 
tect the telescope a "wind watchdog" program stows 
the telescope pointing at the zenith when winds exceed 
steady 8.9 m • s _1 (20 mph) or gusts above 13.4 m • s _1 
(30 mph). The telescope remains stowed until the wind 
speed has remained below these thresholds for 1 hour. 

Observations at zenith angles < 20° are discarded be- 
cause the telescope is unable to track fast enough in 
azimuth to match the sidereal rate near zenith. The 
scheduling algorithm avoids scheduling sources for ob- 
servation at these zenith angles, so few observations are 
lost. Observations at solar or lunar elongations less than 
10° are also discarded. The scheduler does not avoid 
these areas of the sky so a small number of observations 
are lost. 

5.1.2. Pointing and Calibration Failures 

An observation is rejected if a pointing offset was not 
obtained within the prior 4800 s, or if the pointing off- 
set measurement immediately preceding the observation 
failed. Occasional scheduling errors resulted in obser- 
vations without adequately measured pointing offsets. 
These observations are discarded. 

An observation is rejected if fewer than two reliable 
calibration procedures using the CAL diode were suc- 
cessfully executed within a two-hour interval centered on 
the time of the observation, or if the difference between 
the largest and smallest CAL diode measurement within 
that interval differ by more than 10%. 

5.1.3. Saturation or Total Power Anomalies 

The total power varies depending on the attenuator 
setting, receiver gain fluctuations, atmospheric condi- 
tions, and the observed zenith angle. Observations that 
indicate saturation or other total power anomalies are 
rejected. Heavy cloud cover or precipitation often causes 
large fluctuations in total power. Such periods are iden- 
tified by inspection of the total power time series and 
manually discarded. Negative flux density measurements 
are indicated by the 95% upper limits on these values. 

5.1.4. Measured Uncertainty 

We reject flux density measurements with anomalously 
large measured uncertainties, u\5 (equation (J5J) ) . How- 
ever, a straightforward cut at a fixed value or a fixed 
multiple of the expected thermal uncertainty introduces 
a bias against larger flux densities. This occurs because 
there are contributions to the measured error that are 
proportional to the flux density of the target radio source, 
such as telescope tracking errors. We therefore apply a 
flux density-dependent threshold and discard flux density 
measurements for which 

<T16 > C\/l + (p-Sl 5 ) 2 . (9) 



The optimal values of £ = 0.0208 and p = 0.2 were es- 
timated from the data to eliminate as many unreliable 
measurements as possible while minimizing flux density 
bias. About 2% of the data is eliminated by this filter. 

5.1.5. Switched Difference 

We also use the switched difference p, defined by equa- 
tion ([7|. to determine whether flux density measurements 
might be contaminated by systematic errors. The ex- 
pected value of p is 0, provided that the ground spillover 
and atmospheric noise in the ant and ref beams are 
identical. Pointing and tracking errors again give flux 
density-dependent contributions to p, so to avoid bias 
against brighter radio sources we flag points where 

p . o (MQ + Ps-Sis) (1Q) 

Again the optimum values of the parameters (/3 = 5, 
pa = 1.148, p s = 0.0682, and p t = 0.0243) are deter- 
mined from the data. This procedure gives consistent 
results across calibration epochs and it discards about 
2% of flux densities, with comparable fractions dropped 
from each epoch. 

5.2. Relative Calibration 

To correct for slow gain fluctuations of the receiver, we 
first divide each flux density measurement by a calibra- 
tion factor measured using the small noise diode cat A 
measurement of the strength of the cal diode is made af- 
ter each pointing observation, and no less than once per 
hour. Because gain fluctuations are slow, the calibration 
factor is averaged over a two-hour window, centered on 
the time of the flux density measurement. If there are 
fewer than two good measurements of the strength of the 
cal diode in that window then the flux density observa- 
tion is discarded. 

Due to gravitational deformation of the telescope 
structure, the antenna gain varies substantially with 
zenith angle. We model this variation with a polyno- 
mial gain curve and scale flux density measurements to 
remove the effect. Additionally, the optimal axial focus 
position varies with zenith angle, as well as solar zenith 
angle and elongation. During observations, the focus po- 
sition is set using a polynomial model of the zenith an- 
gle variation and a correction is applied during calibra- 
tion using a more complete model that accounts for solar 
zenith angle and elongation. 

The combined effect of these corrections is a factor, 
K re i, that is computed for each flux density measurement. 

5.3. Absolute Calibration 

We divide our observation period into epochs char- 
acterized by a consistent ratio between the calibration 
diode and feed horn inputs to the receiver. This ra- 
tio might change if, for example, the signal path is dis- 
connected and reconnected for maintenance, resulting in 
a slight change in loss along one path. Within a sin- 
gle epoch, the ratio of the calibration diode signal to a 
stable astronomical source should therefore be constant. 
Table [2] lists the epochs we have used in our analysis. 
Absolute calibration is applied to each epoch separately. 
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TABLE 2 
Absolute calibration epochs. 



Epoch 


MJD 


Date 


MJD 


Date 


1 


54466 


2008 Jan 01 


54753 


2008 Oct 14 


2 


54753 


2008 Oct 14 


54763 


2008 Oct 24 


3 


54763 


2008 Oct 24 


55197 


2010 Jan 01 



For each epoch, a calibration factor is determined 
from regular observations of the primary calibrator, 
3C 2 86. We adopt the spectral model and coefficients 
from iBaars et al.l (|1977t) . At our 15 GHz center fre- 
quency, this yields 3.44 Jy, with a quoted absolute uncer- 
tainty of about 5%. The calibration factor for epoch i, Atj, 
is the ratio of the adopted flux density for the calibrator 
to the weighted mean of the observations: 



(LSis-o^ 2 )/^^ 2 )' 

where S[ 5 and er 15 denote the flux densities for the cali- 
brator with only the relative calibration applied. 

The total calibration factor for a flux density mea- 
surement in equation Q is then and 
reflects both relative and absolute calibration. Cross- 
checks of our calibration against 14.6 GHz observations 
of a number of common sources observed with the Ef- 
felsberg 100 m telescope through the F-GAMMA project 
confirm the overall accuracy of our flux density scale. 

5.4. Uncertainties in Individual Flux Density 
Measurements 

In a perfect observing system with no sources of sys- 
tematic error the uncertainties in the flux density mea- 
surements would be given by the thermal noise on each 
observation. In practice there are many sources of sys- 
tematic error, including the effects of weather and the 
atmosphere, mis-pointing due to wind, and focus errors. 
Many of these are correctly identified and accounted for 
in the automatic and manual editing described in § 15.11 
However, even after flux density measurements affected 
by these problems are filtered out there remain many ob- 
servations that are significantly affected by systematic er- 
rors. Such systematic errors can lead to significant errors 
in the measurement that are not reflected in the thermal 
noise of the observation and can give rise to bad flux den- 
sity measurements with small thermal errors. This leads 
to "outliers" in the light curves, i.e. points which do not 
lie close to the level determined from interpolation of ad- 
jacent observations and which have small errors. The 
task of identifying and eliminating or allowing for the 
wide variety of systematic errors leading to such outliers 
is challenging and time-consuming. Great care must be 
taken not to assume that the behaviour of the source 
is known, and hence to eliminate a real and potentially 
extremely interesting flux density variation. 

We first apply an error model to determine the uncer- 
tainty of each flux density measurement: 

^tota! = *? B + (e ■ S 15 f + (n ■ V) 2 , (12) 

w hich is an extensio n of the model described 
in lAngelakis et al.l (|2009f) . The first term represents the 



TABLE 3 
Error Model Parameter Values 



Parameter 


Pointing Calibrator 


Early a 


Late a 


e 


0.0057 


0.0200 


0.0135 


1 


3.173 


3.173 


3.173 



a The "early" model applies prior to MJD 54906 (2009 
March 16). 



measured scatter during the flux density measurement. 
This includes thermal noise, rapid atmospheric fluctua- 
tions, and other random errors. The second term adds an 
uncertainty proportional to the flux density of the source. 
This term allows for pointing and tracking errors, varia- 
tions in atmospheric opacity, and other effects that have 
a multiplicative effect on the measured flux density. In 
the third term tf> is the switched power, defined by equa- 
tion (J6j) . This term takes account of systematic effects 
that cause the A-B segment of the flux density measure- 
ment to differ from the C-D segment, such as a pointing 
offset between the A and D segments, or some rapidly 
varying weather conditions. 

The error model is defined by the two parameters, e 
and n, whose values must be determined from the obser- 
vations. Because e describes the error contribution due 
to pointing errors, its value depends on whether a source 
is used as a pointing calibrator. Furthermore, for non- 
pointing sources, e is found to differ between the schedul- 
ing algorithms used before and after MJD 54906 (2009 
March 16). The parameter, r), is found to be adequately 
described by a single value for all sources and all epochs. 
The adopted values are given in Table [3] 

For pointing sources, both e and 7/ were estimated si- 
multaneously using the stable flux calibrators 3C 286, 
3C 48, 3C 161, and DR 21. Due to systematic errors, 
these sources and other stable-flux density calibrators 
show long-term variations of 1-2% so we fitted a 7th- 
order polynomial to remove this trend from each source, 
then computed the residual standard deviation, median 
flux density, the rms, and mean ip for each source, then 
used these to fit the error model parameters. 

To determine the error model parameter e for ordi- 
nary sources, we selected 100 sources that exhibited little 
variation or slow, low-amplitude variations in flux den- 
sity, between the start of our program and MJD 55048 
(2009 August 5). This interval was split into two periods, 
"early" and "late," at MJD 54906 and this procedure 
was separately applied to each period. For each light 
curve, we fitted and removed a second-order polynomial 
trend, then iteratively removed outlier data points with 
residuals greater than three standard deviations. We re- 
peated the fitting and outlier removal until no further 
outliers were removed and we discarded any source with 
fewer than 10 remaining data points (retaining 94 and 88 
sources in the early and late periods, respectively). From 
the surviving points in each light curve, we computed the 
median and the rms flux densities, and the standard de- 
viation of the residuals. We then fitted equation (|12l) to 
these data, omitting the r\ term. The data and the error 
model results for the early period are shown in Figure |31 
We then adopted the same value of rj for these sources 
as was determined for pointing sources. 
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Fig. 3. — Residual standard deviation (gray points) and fitted e- 
only error model values (black crosses) for ordinary sources in the 
early (MJD < 54906) period. The fit in the late period is similar. 
A single high-flux density data point was omitted to limit the scale. 




100 200 300 400 500 600 
Days since MJD 54466 (2008 Jan 1) 

Fig. 4.— Normalized flux densities for 3C 274 (top), DR 21 
(center), and 3C 286 (bottom) after outlier removal. Each light 
curve is normalized by its median. The black line in each plot is 
the spline fit to the combined data. 



5.4.1. Long-Term Trends in 3C 286, 3C 27^, and DR 21 

After carrying out the above editing and calibration 
steps we returned to the residual 1-2% long-term (~ 
6-month) variations in the light curves for stable-flux- 
density calibration sources. We chose 3C 286, 3C 274, 
and DR 21 for this study because they are well-known 
to be stable on timescales of many years. The fractional 
variations in flux density of these objects are shown in 
Figure 2] and are clearly correlated, indicating the pres- 
ence of an unidentified source of multiplicative system- 
atic error. For each of these sources, we removed 2-a 
outliers in a 100-day sliding window and normalized the 
resulting data by the median flux density. We then com- 
bined the data for all three sources and fitted a cubic 
spline to the result. 

We apply the corresponding correction to all light 
curves in our program by dividing each flux density by 
the value of this spline. Figure [5] shows the residuals 
for the three fitted sources after dividing out the spline 
fit. The 1% residual variation that remains is the level of 
systematic uncertainty after correction for this long-term 
trend. 

5.5. Scaling of the Non-Thermal Error 
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Fig. 5.— Normalized flux densities for 3C 274 (top), DR 21 
(center), and 3C 286 (bottom) after dividing by the spline fit to 
remove long-term systematic trends. 

The reported error for each flux measurement has two 
qualitatively different components as described in section 
15.41 The first component is directly obtained during the 
flux measurement and it represents random errors such as 
thermal noise and rapid atmospheric fluctuations, while 
the second is introduced to take into account other, flux- 
density-dependent effects. This error model requires the 
determination of constant factors, which we have called 
e and rj, and which have been assumed to be source inde- 
pendent. However, many sources exhibit coherent long- 
term variations with scatter about those clearly smaller 
than what would be expected as a result of the quoted 
errors. This is a direct indication that in certain cases 
the simple assumption of source-independent e and r] re- 
sulted in overestimated errors. 

To correct these constant scale factors on a source-by- 
source basis, we have used cubic spline fits and required 
a x 2 P er degree of freedom to be one for the residuals. 
Due to the large number of sources and the requirement 
of an uniform and consistent method for all the sources, 
an automatic method was developed for this procedure. 
For each source we can in principle use a range of num- 
ber of polynomial sections to construct a spline fit. We 
construct a spline fit for each possible number of poly- 
nomial sectionsQ An outlier rejection filter which uses 
a cubic spline fit with a small number of knots is used 
to fit the light curve. Points with absolute residuals in 
the largest 5% are not used for the following stage of the 
fitting procedure. Not all the fits are acceptable, as some 
cases will have correlated residuals or a large departure 
from normality. Acceptable fits are selected by using two 
statistical tests: Lilliefors test for normality and the runs 
test for randomness^ Only the fits for which both null 
hypotheses cannot be rejected at the 10~ 3 level are con- 
sidered acceptable. For each acceptable fit, a scale factor 
that makes the \ 2 P er degree of freedom equal to one is 
calculated. Among the scale factors for all the acceptable 
fits, the median scale factor is selected as the final correc- 
tion. The value of the scale factor is not very sensitive 
to the exact number of polynomial sections. A typical 
example of the behavior of the scale factor is shown in 

1 We use the MATLAB Spline Toolbox function spap2, which 
automatically selects the positions of the knots for the spline. 

2 We have used the implementation of both tests that are part 
of the MATLAB Statistics Toolbox. 
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Fig. 6. — Example of the error scale factor correction using data 
for J0046+3900. The two upper panels show the light curve with 
the original (left) and corrected (right) error bars (gray points) 
and a typical spline fit (black line). The bottom left panel shows 
the residuals from the spline fit using the corrected error bars. 
In the bottom right panel, the \ 2 P er degrees of freedom (solid 
gray line) and correction factor (solid black line) are shown, with 
black circles marking the correction factors for fits that pass the 
acceptance tests, and a dashed line showing the adopted correction 
factor for the source. 



Figure [6l 

We have thus only rescaled the non-thermal part of 
the errors (the Si 5 and ip terms in equation (|12p). and 
only for those sources for which the resulting correction 
factor was smaller than one (i.e. the rescaling would 
result in smaller errors) . The latter choice was made for 
two reasons. First, a correction factor larger than one 
simply indicates that the spline fit cannot provide an 
adequate description of the data. This may result from 
a light curve more variable than can be fit by spline with 
a given number of knots, so such a correction could mask 
real variability. Only the reverse is cause for concern — 
when the spline fit is too good a fit, given the quoted 
errors. Second, this choice ensures a smooth transition 
between scaled and non-scaled errors, as the transition 
point (correction factor equal to one) is equivalent to no 
error scaling. 

6. RESULTS AND ANALYSIS 

6.1. Monitoring Program Statistics 

Our target cadence was two flux density measurements 
per source per week, or about 200 measurements per 
source in the first two years of the program. Our re- 
sulting average effective cadence for CGRaBS sources is 
about 134 measurements per source in the first two years 
of the program. The efficiency compared to our nominal 
cadence is 67%. 

6.2. Light Curves 

Light curves for the CGRaBS program sources are 
shown in Figures 7.1-7.1158. Table [S] lists the filtered 
and calibrated 15 GHz flux density measurements that 
result from the procedure described above. Regular up- 
dates to the data set, including data for sources outside 
the core sample released in this paper, are available from 



the program website^ 

6.3. Source Variability 

In this section we discuss the variability amplitude 
observed in each source in our sample. The variabil- 
ity properties of our sources in the time and frequency 
domains, as quantified by measures such as the power 
spectrum and autocorrelation function, and correlation 
with gamma-ray data to identify and measure time lags, 
will be discussed in detail in a forthcoming publica- 
tion (|Max-Moerbeck et al.ll2011[ ). 

The questions of the variability amplitude of a source 
and the confidence with which this can be measured are 
complex ones and have been traditionally addressed us- 
ing a variety of measures and test s, such as the vari- 
ability in dex (e.g.. lAller et al.lll992T) : the fluctuation in- 
dex (e.g..lAller et al.ll2003[ ): the modulation index (e.g., 
iKraus et al.l I2003D; the fractional variability amplitude 
(e.g.. lEdelson et al.H2002T; ISoldi et al JI2008Q : and x 2 tests 
of a null hypothesis of non- variability. Each of these tools 
provides different insights to the variability properties of 
sources and is sensitive to different uncertainties, biases, 
and systematic errors. For example, the variability in- 
dex, defined as the peak-to-trough amplitude change of 
the flux, is a measure of the amplitude of the variability 
of a source: 



V = 



(S n 



(S n 



(S n 



(S n 



(13) 



where S max and S m - ln are the highest and lowest mea- 
sured flux densities, respectively, and a max and <J m i n are 
the uncertainties in these measurements. Although the 
definition is constructed to account for the effect of mea- 
surement uncertainties, the quantity is well-defined only 
when variability is significantly greater than measure- 
ment errors, and it can yield negative values for sources 
with low signal-to-noise ratios or little intrinsic variabil- 
ity. In addition, it is very sensitive to outliers and is 
not robust against random Gaussian excursions from the 
mean. Such excursions are to be expected for sources 
that are regularly monitored over long periods of time: 
even non-variable sources are likely to have at least one 
pair of 2-cr high and low measurements after being ob- 
served more than 100 times, as is the case for most 
sources in our sample. 

An associated measure of variability amplitude is the 
modulation index, defined as the standard deviation of 
the flux density measurements in units of the mean mea- 
sured flux density, 



?™data — 



1 v-^ ( q L 



N S 



(14) 



The modulation index has the advantage that it is always 
non- negative and more robust against outliers. However, 
it still represents a convolution of intrinsic source varia- 
tion and observational uncertainties: a large modulation 
index could be indicative of either a strongly variable 
source or a faint source with high uncertainties in indi- 
vidual flux density measurements. For this reason, the 

3 http:/ /www. astro. caltech.edu/ovroblazars 
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TABLE 4 

Program Source Data 3 - 



Name 






_b 
z 


Optical 
Class b 


IN U1T1 

Obs 


flag 


■So 
(Iv) 


m c 
<%) 


3C 48 


01 h 37 m 41?30 


33°09'35('l 






281 





1.805 ± 0.002 


1.8 ±0.1 


3C 161 


U6 27 llnlz 


— 05 53 05. 2 






191 





2.081 ± 0.005 


3.3 ± 0.2 


3C 274 


12 h 30 m 49f42 


12°23'28('0 






252 


5 


26.329 ± 0.021 


0.9 ±0.1 


3C 286 


13 h 31 m 08?29 


30°30'33?0 






232 


5 


3.438 ± 0.003 


0.7 ± 0.1 


DR 21 


20 h 39 m 01?20 


42°19'32'9 






282 


5 


19.024 ± 0.011 


0.7 ± 0.1 


J0001+1914 


00 h 01 m 08?62 


19°14'33'.'8 


3.100 


FSRQ 


165 





0.282 ± 0.003 


11 5+ ' 8 


J0001-1551 


00 h 01 m 05?33 


-15°51'07'.'l 


2.044 


FSRQ 


159 





0.213 ± 0.002 


8 9 +0 - 7 


J0003+2129 


00 h 03 m 19?35 


21°29'44'.'4 


0.450 


AGN 


169 





0.087 ± 0.001 


7 n+0.8 
' -0.7 


J0004+2019 


00 h 04 m 35?76 


20°19'42'.'2 


0.677 


BLL 


175 





0.327 ±0.003 


12.5lg' 8 7 


J0004+4615 


00 h 04 m 16?13 


46°15'18'/0 


1.810 


FSRQ 


154 





0.181 ± 0.006 


38.5±S;| 


J0004-1148 


00 h 04 m 04=92 


-11°48'58'.'4 




BLL 


106 





0.720 ± 0.011 


15.511;? 



a Only the first few rows are shown here; the complete version of this table is available in the electronic 
version of the journal. 

b For the CGRaBS sample, redshift and optical classifications are repeated here from lHealev et al. (2008) for 
convenience. 

c The number of observations that survived data editing and were used in our variability analysis. 

d Variability analysis flag. A value of indicates a non-zero intrinsic modulation is found; 1 indicates the 

source is non-variable; 2 indicates insufficient observations for variability analysis; 3 indicates flux density too 

faint for variability analysis; 5 indicates the source was a calibrator used in the spline fit to remove long-term 

trends. 

e Quoted errors are 1-a uncertainties. Values for non-variable sources indicate 3-cr upper limits and no 
uncertainties are quoted for m or Sq. 



TABLE 5 

15 GHz Flux Densities 3 



Source 


MJD 


Flux Density (Jy) 


J0001- 


1551 


54471.051377 


0.244 ± 0.008 


J0001- 


1551 


54474.042836 


0.232 ± 0.007 


J0001- 


1551 


54478.032303 


0.221 ± 0.008 


J0001- 


1551 


54480.026840 


0.238 ± 0.011 


J0001- 


1551 


54484.015903 


0.229 ± 0.008 



a Only the first few rows are shown here; the com- 
plete version of this table is available in the elec- 
tronic version of the journal. 



correct interpretation of results on the modulation index 
requires that measurement errors and the uncertainty in 
m due to the finite number of flux density measurements 
be properly accounted for. 

One method that has been widely used to evaluate the 
information encoded in variability measures is to evaluate 
each measure for a set of constant-flux-density calibra- 
tors, which are known to have a flux density constant in 
time and have been observed with the same instrument 
over the same periods of time. The value of the vari- 
ability measure obtained for the calibrators is then used 
as a threshold value, so that any source with variability 
measure equal to or lower than that of the calibrators 
is considered consistent with being non-variable. How- 
ever, a variability measure value higher than that of the 
calibrators is a necessary but not sufficient condition for 
establishing variability: calibrators are generally bright 
sources, with relative flux density measurement uncer- 
tainties typically lower than the majority of monitored 
sources; additionally, variability measures are affected by 
the sampling frequency, which is not necessarily the same 
for all monitored sources and the calibrators. 



Alternatively, the significance of variability in a given 
source can be established through tests (such as a x 2 
test) evaluating the consistency of the obtained set of 
measurements with the hypothesis that the source was 
constant over the observation interval. However, such 
tests provide very little information on sources for which 
statistically significant variability cannot be established, 
as they cannot distinguish between intrinsically non- 
variable sources and sources that could be intrinsically 
variable but inadequately observed for their variability 
to be revealed. 

Here we propose a new index for characterizing source 
variability: the intrinsic modulation index m, which is 
the intrinsic standard deviation of the distribution of 
source flux densities in time, o~q, measured in units of 
the intrinsic source mean flux density, Sq . Here the term 
"intrinsic" is used to denote flux densities and variations 
as would be observed with perfectly uniform sampling of 
adequate cadence and zero observational error: 



In this way, m is a measure of the true amplitude of 
variations in the source, rather than a convolution of 
true variability, observational uncertainties, and effects 
of finite sampling. Observational uncertainties and finite 
sampling will, of course, affect the accuracy with which m 
can be measured. The purpose of the analysis described 
in this section is to derive a best estimate of m, as well 
as an estimate of the uncertainty in our measurement of 
this quantity. For sources with m within 3c from zero, 
the 3<t upper limit on m will be evaluated. 

6.3.1. A Likelihood Analysis to Obtain the Intrinsic 
Modulation Index 
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Fig. 7. — 15 GHz light curves for calibrators and CGRaBS program sources. Light curves for all sources are available in Figures 7.1—7.291 
in the electronic version of the journal. 



For the purposes of our analysis, we will assume that 
the "true" flux densities for each AGN are normally dis- 
tributed, with mean So, standard deviation uq, and in- 
trinsic modulation index m = ao/So- We have N mea- 
surements of the flux density, Sj, each of which has an 
associated observational uncertainty, also assumed Gaus- 
sian, tjj. 

Let us assume that at a moment of observation, a 
source has a "true" flux density St- The probability 
density to observe a value near Sj if the observational 
uncertainty is tjj is 



1 



'2tt 



exp 



2a) 



(16) 



In addition, the probability density that the true source 
flux density at one of the moments of observation is near 
Si if the source flux densities are distributed normally 
with mean So and standard deviation ao is 



p(S t ,S ,a ) = 



1 



exp 



(S t - S f 



2al 



(17) 



Therefore, the likelihood of observing one flux density Sj 



with uncertainty Oj from the particular source is 



= / dS t 

all S t 



exp 



exp 



(18) 



which amounts to calculating the probability to observe 
Sj through any possible true flux density value St ■ If the 
limits of integration above are taken to be from St = — oo 
to S t = oo then the integral ha s an analytic solution (see, 
e.g. JVenters fc Pavlidoull2007D : 



: exp 



(S 3 - Sq) 2 
2K + cr 2 ) 



(19) 



The likelihood for N observations (Sj,<Tj) for j 
l,..JVis 



C(So, ao) 




N 



1 



exp 



1 N 



3=1 



{S 3 -Spf 
a] + al 



(20) 
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Fig. 8. — 1, 2, and 3<r contours of the joint likelihood £(So> m ) 
for blazar J1243-0218. 



The intrinsic standard deviation ao can be eliminated in 
favor of the intrinsic modulation index, 



<jq = mS , 



(21) 



so that 



C(S ,m) 




2tt(to 2 S 2 + a 2 ) 



exp 



1 (Sj - S qY 

9. ^—^ rr 2 -J- rrf- 



2 — a- 

2 = 1 3 



o2 

o 



(22) 



This likelihood is symmetric about to = 0, as to only 
enters through its square. For this reason, this formalism 
can guarantee non-negative intrinsic modulation indices 
without loss of information. 

Maximizing the joint likelihood C(So,rn), we can de- 
rive maximum-likelihood estimates of So and to. Isolike- 
lihood contours containing 68.26%, 95.45%, and 99.73% 
of the total volume under the joint likelihood surface de- 
fine the 1, 2, and 3a contours respectively (see Figure |5] 
for an example in the case of J1243— 0218, whose light 
curve is shown in Figure [7]). The maximum- likelihood 
Gaussian for the distribution of flux densities for the 
same object is compared to the histogram of measure- 
ments in Figure |9l Note that the maximum-likelihood 
Gaussian is narrower than the histogram; this behav- 
ior is expected, as the histogram is a representation of 
measurements sampling the underlying distribution with 
finite error. The typical magnitude of the latter for the 
particular source is shown in Figure |9] with the blue ar- 
rows, and it is indeed comparable with the difference 
in width between the maximum-likelihood Gaussian and 
the histogram. 

To derive the most likely value of fn and the associated 
uncertainties regardless of the true value of So, we inte- 
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Fig. 9. — Maximum-likelihood Gaussian model for the flux den- 
sity distribution (dashed line) . plotted over the histogram of mea- 
sured flux densities (solid line) for blazar J1243— 0218. The arrow 
indicates the size of the typical measurement uncertainty. 



grate So out of £(So,m), and obtain the marginalized 
likelihood as a function of only to: 



L(m) — 



dSo So 

all So 



exp 




N 



-y 



(Sj - So) 2 



i 2 si 



(23) 



Then, the value of fn that maximizes the marginalized 
likelihood is our best estimate of it, and the 1-a uncer- 
tainty on the modulation index can be found by locating 
the isolikelihood to— values toi and TO2 for which 



C(m\) = C{m2) 



and 



r 2 jfc/jrAffiYi 

Jo°° £(m)dm 



= 0.6826. 



(24) 



(25) 



Note that the upper and lower la errors are not gener- 
ally symmetric in our formalism. Similarly, 2a and 3a 
ranges can be derived by substituting the right-hand- 
side of equation ((25|) by 0.9545 and 0.9973, respectively. 
The marginalized likelihood, best-estimate to, and the 
la and 2a m ranges for blazar J1243— 0218 are shown in 
Figure [KB 

If the maximum- likelihood to is less than 3a away from 
to = 0, we consider that statistically significant variabil- 
ity cannot be established. In these cases, we calculate 
the 3a upper limit on to, which is defined as the value 
to 3 for which 



J 3 C(m)drfi 
J °° £(m)dm 



0.9973. 



(26) 
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Fig. 10.— Marginalized likelihood £(m) for J1243-0218 (solid 
curve). Dashed vertical line: best-estimate m; dotted vertical lines: 
la m range; solid vertical lines: 2<r m range. 

The use of the intrinsic modulation index and the like- 
lihood analysis we have employed to estimate it have the 
advantage of offering a way to obtain information about 
the intrinsic variability of the source, deconvolved from 
observational errors in individual flux density measure- 
ments and the effects of finite sampling, while providing 
strictly defined 1, 2 and 3cr uncertainties for our estimate 
of to (essential when conducting population studies), and 
upper limits for to when variability cannot be established 
at a > 3<7 confidence. However, our choice carries its own 
caveats. 

(1) Model- dependence of the likelihood analysis. A 
functional form has to be assumed for the intrinsic dis- 
tribution of flux densities (here we have assumed it to 
be Gaussian), resulting in a loss of generality. The va- 
lidity of this assumption can be tested by comparing 
the maximum-likelihood intrinsic flux density distribu- 
tion to the histogram of measured flux densities, to eval- 
uate whether the maximum-likelihood flux density dis- 
tribution is a reasonable description of the data (mod- 
ulo observational uncertainties). This is indeed the case 
for many, although not all, of our sources. An example 
of a source well described by the maximum-likelihood 
flux density distribution is shown in Figure |H1 Other 
sources however show bimodality, and the distribution 
of measured flux densities in these cases could be better 
described by, for example, a double Gaussian. An ex- 
tended likelihood analysis that does explicitly account for 
bimodality and calculates not only the intrinsic modula- 
tion index but also duty cycles and flaring-to-quiescent 
flux density ratios will be presented in an upcoming pub- 
lication. For the purposes of this work we have con- 
firmed that, even when a single Gaussian is not an ad- 
equate description of the flux density distribution, the 
intrinsic modulation to index is well-correlated, within 
uncertainties, with the modulation index TOdata derived 
from the data (equation (fT4)) ). which, although con- 
taminated by observational uncertainties, is completely 



model- independent (see next section). 

(2) Assumption of unbiased sampling. Our analysis as- 
sumed that the flux density values we have not sampled 
are not correlated with each other. This assumption is 
poor in the case of lengthy outages, as well as for in- 
creased cadence for any single epoch. In our analysis we 
have disregarded the additional data taken during epochs 
of increased cadence for specific objects (during, for ex- 
ample, campaigns to constrain intra-day variability). 

(3) Leakage of probability density to negative flux den- 
sities. In certain cases, extending the integration over 
intrinsic flux densities from — oo to oo (in equation (|18jl) 
to simplify the mathematical manipulations leads to un- 
acceptable leakage of probability density to unphysical 
domain of negative true flux densities. This approxi- 
mation is adopted for numerical efficiency, since in this 
case the likelihood can be expressed analytically with- 
out the need to perform multi-dimensional integrals for 
every object in our large sample. For most objects in 
our sample the leakage to the negative flux density do- 
main is negligible. The error introduced in this way be- 
comes important only for very dim AGN (because the 
peak of the St distribution is very close to zero) or very 
variable AGN (because of very long tails in the St dis- 
tribution). None of the CGRaBS sources in our sample 
are dim enough for the first effect to be a problem, and 
very few are variable enough: for to ~ 0.5, about 3% of 
the "true flux density" probability density leaks to neg- 
ative values, with the problem becoming more severe for 
more variable sources; only four CGRaBS sources have 
to > 0.5. 

6.3.2. Variability Analysis — Results 

In Figure [TT] we plot the intrinsic modulation index 
to and associated 1-cr uncertainty against the intrin- 
sic, maximum-likelihood average flux density, So, for all 
sources in the OVRO monitoring program (CGRaBS and 
non-CGRaBS sources). The error bar on Sq corresponds 
to the 1-tr uncertainty in mean flux density, calculated 
from the joint likelihood (equation (|22j) ) marginalized 
over to. CGRaBS sources are shown in black, while non- 
CGRaBS sources are shown in red. 

Variability could only be established at the 3-cr con- 
fidence level or higher for 1146 out of 1158 CGRaBS 
blazars in our sample. For this study, we considered only 
sources for which at least 3 flux densities were measured, 
a positive mean flux density > 2a from zero was found, 
and at least 90% of the individual flux density measure- 
ments were > 2a from zero. These criteria excluded one 
source. For the other 11 sources we have calculated 3-a 
upper limits for to. We plot these upper limits with blue 
triangles. We also plot, with magenta triangles, 3-er up- 
per limits for non-CGRaBS sources that were consistent 
with non- variable and had enough flux density measure- 
ments. 

Calibration sources 3C 286, DR 21, and 3C 274 are 
shown in green. Although these sources are the least vari- 
able (as expected) of all sources in which variability can 
be established and a non-zero to can be measured, to for 
these sources is finite and measurable. This means that 
some residual long-term variability remains in our cali- 
brators beyond what can be justified by statistical errors 
alone. This could conceivably result from true calibra- 
tor source variation, but more likely reflects incomplete 
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Fig. 11. — Intrinsic modulation index m and associated l-cr 
uncertainty, plotted against intrinsic maximum-likelihood average 
flux density, So, for all sources in the OVRO monitoring program 
which have enough (more than 3) acceptable, non-negative flux 
density measurements. Black points: CGRaBS sources; red points: 
non-CGRaBS sources; green points: calibrators; blue triangles: 3- 
a upper limits for CGRaBS sources for which variability could not 
be established at > 3cr confidence level; magenta triangles: as blue 
triangles, for non-CGRaBS sources. The error bar on So corre- 
sponds to the 1-tr uncertainty i n m ean flux density, calculated from 
the joint likelihood (equation J22)) marginalized over m. Variable 
CGRaBS sources outside the yellow and cyan shaded areas are used 
in the population studies of § 16.3.31 

removal of small-amplitude calibration trends. Because 
TO < 1% for these three sources, we quote a systematic 
uncertainty Am syst = 0.01 for the values of the intrinsic 
modulation index we produce through our analysis. 

To ensure that our population studies are not affected 
by this residual systematic variability, in all analyses dis- 
cussed in § 16.3.31 only sources with m > 0.02 will be 
used, so that we remain comfortably above this 1% sys- 
tematic uncertainty limit. In addition, for sources with 
So < 60 mjy, the number of sources for which variability 
can be established is of the same order as the number of 
sources (both CGRaBS and non-CGRaBS) for which we 
could only measure an upper limit, and these upper lim- 
its are very weak and non-constraining. For this reason, 
we also exclude from our population studies of § 16.3.31 all 
sources below So = 60 mjy. The part of the parameter 
space excluded due to these two criteria is shown in Fig- 
ure [TT] as the yellow shaded area bounded by the solid 
black lines. 

For 5*o > 0.4 Jy, no obvious correlation between 
flux density and modulation index is apparent, and no 
CGRaBS sources exist with upper limits above our cut 
of m = 2%. However, for sources with So < 0.4 Jy, 
there is an absence of points in the lower-left corner of 
the allowed parameter space defined by the thick solid 
lines: for faint sources, we can only confidently estab- 
lish variability if that variability is strong enough. The 
effect disappears for variability amplitudes greater than 
about 6%. In addition, there are no CGRaBS upper lim- 
its higher than 6% for sources brighter than 60 mjy. We 
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Fig. 12. — Intrinsic modulation index m and associated l-cr 
uncertainty, plot ted against the "raw" modulation index, m<j ata , 
of equation J 1 41) as black points with brown error bars. The 
m = nidata line is shown in blue. Green triangles are the 3-cr up- 
per limits of sources for which variability could not be established. 
Calibrators are plotted in red. 



conclude that we are able to measure variability at the 
level of 6% or higher for any CGRaBS source brighter 
than 60 mjy. 

To ensure that our population studies are not affected 
by our decreased efficiency in measuring variability in 
sources with 60 mjy < S < 0.4 Jy and 2% < m < 6%, 
we will also exclude this part of the (So,m) parameter 
space from our analysis in § 16.3.31 The part of the pa- 
rameter space excluded due to these criteria is shown in 
Figure [Tl] as the cyan shaded area. 

The single point with m > 1.0 near the upper left cor- 
ner of the plot is Cygnus X-3, a galactic microquasar not 
part of our core program, which indeed exhibits a strong 
flare during which the flux density increases by about an 
order-of-magnitude over its average in a short period of 
time. Although in this case our underlying assumption 
of the flux density measurements being consistent with 
a single Gaussian distribution is clearly not a sufficient 
description of the data (rather, this case is much better 
described by a bimodal distribution), the qualitative re- 
sult of variability with very high amplitude in this source 
is robust. 

In Figure [12] we plot the intrinsic modulation index m 
and associated l-cr uncertainty against the "raw" modu- 
lation index mdata of equation (fTl)) . The m = mdata line 
is shown in blue. Green triangles are the 3-cr upper limits 
of sources for which variability could not be established. 
Calibrators 3C 286, DR 21, and 3C 274 are plotted in 
red. Since apparent variability due to the finite accuracy 
with which individual flux densities can be measured has 
been corrected out of m, the expectation is that devia- 
tions from the m = mdata will be more pronounced for 
sources that are not intrinsically very variable (so that 
the scatter in the flux density measurements is appre- 
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Fig. 13. — Histogram of maximum-likelihood intrinsic modula- 
tion indices m, for the 456 CGRaBS blazars with So > 400 mjy. 
The dashed line represents an exponential distribution with (m) = 
0.091. 



ciably affected, and even dominated, by measurement 
error). In addition, deviations are expected to be be- 
low the line, as m should be smaller than mdata- Both 
these expectations are verified by Figure [TTJ Note that 
upper limits need not satisfy this criterion, as the "true" 
value of the modulation index can take any value below 
the limit. Upper limits above the blue line are weak, 
indicating that the reason variability could not be estab- 
lished is the poor sampling or quality of the data, and 
not necessarily a low intrinsic variation in the source flux 
density. 

For the 456 CGRaBS objects which have S > 400 mjy 
and for which variability can be established, we plot, in 
Figure 1131 a histogram of their intrinsic modulation in- 
dices m normalized so that the vertical axis has units of 
probability density. The dashed line represents an expo- 
nential distribution of mean (m) = 0.091 which, as we 
can see, is an excellent description of the data. Moti- 
vated by this plot, we will be using the monoparametric 
exponential family of distributions: 



f(m)dm = — exp 
m 



m 
m 



dm 



(27) 



with mean mo and variance m§, to characterize various 
sub-samples of our blazar sample. 

6.3.3. Variability Analysis — Population 
Studies — Formalism 

We now turn our attention to whether the intrinsic 
variability amplitude at 15 GHz, as quantified by m, cor- 
relates with the physical properties of the sources in our 
sample. To this end, we will determine the distribution 
of intrinsic variability indices m for various subsets of 
our monitoring sample, and we will examine whether the 



various subsets are consistent with being drawn from the 
same distribution. 

We will do so using again a likelihood analysis. We will 
assume that the distribution of rn in any subset is an ex- 
ponential distribution of the form given in equation (I27p . 
Since distributions of this family are uniquely described 
by the value of the mean, mo, our aim is to determine 
mo, or rather the probability distribution of possible mo 
values, in any specific subset. 

The likelihood of a single observation of a modulation 
index rrii of Gaussian uncertainty <7i drawn from an ex- 
ponential distribution of mean mo is 




where, to obtain the second expression, we have com- 
pleted the square in the exponent of the integrand. The 
last integral can be calculated analytically, yielding 
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(29) 



If we want (as is the case for our data set) to implement 
data cuts that restrict the values of rrii to be larger than 
some limiting value m; , the likelihood of a single obser- 
vation of a modulation index fS, will be the expression 
above multiplied by a Heaviside step function, and renor- 
malized so that the likelihood ^i, C uts to obtain any value 
of rrii above mi is 1: 



£i,cuts mil 



H(m,j - mi)£j 
f— drriiii 

Jm,—mi L L 



(30) 



This renormalization enforces that there is no probability 
density for observed events "leaking" in the parameter 
space of rejected rfTi values. In this way, it "informs" the 
likelihood that the reason why no objects of rrii < mi 
are observed is not because such objects are not found in 
nature, but rather because we have excluded them "by 
hand." 

The integral in the denominator is analytically calcu- 
lable, 



drriili = - 
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1 +erf 



-1 - erf 
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The likelihood of N observations of this type is 

N 

£(m ) = Y\_£i, cu ts[mi] ■ 



(31) 



(32) 
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If we wish to study two parts of the Sq parameter space 
with different cuts (as in, for example, Figure [Til where 
we have a cut of mi = 0.02 for So > 0.4 Jy, and a differ- 
ent cut of m u = 0.06 for 0.06 Jy < Sq < 0.4 Jy), we can 
implement this in a straight-forward way, by consider- 
ing each segment of the Sq parameter space as a distinct 
"experiment", with its own data cut. If the first "exper- 
iment" involves Ni objects surviving the mi cut, and the 
second "experiment" involves N u objects surviving the 
m u cut, then the overall likelihood will simply be 



N, 



c(m ) = l[e 

i,CUts[^z] I I ^2, CUtS [^li] • 



(33) 



Maximizing equation (|33[) we obtain the maximum- 
likelihood value of mo, »no,maxL- Statistical uncertainties 
on this value can also be obtained in a straight-forward 
way, as equation (I33|) . assuming a flat prior on mo, gives 
the probability density of the mean intrinsic modulation 
index mo of the subset under study. 

6.3.4. Variability Analysis — Population Studies — Results 



Here we apply the formalism introduced in § 16.3.31 to 
examine whether the intrinsic modulation index m cor- 
relates with the physical properties of the sources in our 
sample. We will be testing whether the distributions of 
m-values in subsets of our monitoring sample split ac- 
cording to some source property are consistent with each 
other. To verify that our analysis does not yield spu- 
rious results, we first discuss two test cases where the 
likelihood analysis should not find a difference in the 
variability properties of the different subsets considered. 

The first case tests whether the data cuts discussed 
in § 16.3.21 are implemented correctly in § 16.3.31 To this 
end, we calculate C(mo) for the set of non gamma-ray- 
loud CGRaBS blazars (blazars not found in 1LAC) in 
our monitoring sample with S > 0.4 Jy, in two different 
ways: first, by applying an m cut at m; = 0.02; second, 
by applying an m cut at mi = 0.06 (a much more aggres- 
sive cut than necessary for the particular bright blazar 
population). The increased value of m; in the second 
case should not affect the result other than by reduc- 
ing the number of data points and thus resulting in a 
less constraining likelihood for mo. This is indeed the 
case, as we see in Figure [T4j where we plot the proba- 
bility density of mo for the two subsets. That the two 
distributions are consistent with each other is explicitly 
demonstrated in Figure [T5l where we plot the probability 
density of the difference between the means mo of the two 
subsets (which is formally equal to the cross-correlation 
of their individual distributions). The difference is con- 
sistent with zero within la. 

The second case tests whether a split according to a 
source property without physical meaning and with the 
same value for the cutoff modulation index m; will yield 
probability densities for the mo that are consistent with 
each other. For this reason, we split the population of 
bright (S > 0.4 Jy) CGRaBS blazars in our monitoring 
sample in two subsets in the following way: we divide 
the RA of each source by 1 min. If the remainder of this 
operation is < 30 s, we include this source in the first 
subsample (depicted by a solid line in Figure [T6]) . If the 
remainder is > 30 s we include the source in the second 
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Fig. 14. — Probability density of mo for the subset of bright 
CGRaBS blazars not found in 1LAC, for two values of the cutoff 
for data acceptance: m; = 0.02 (solid line), and m; = 0.06 (dashed 
line). The two distributions are consistent with a single value. 




Fig. 15. — Probability density of the difference between t he m ean 
modulation index mo for the two sets considered in Figure [l4l The 
difference is consistent with zero within lc. 



subsample (depicted by a dashed line in Figure [15]). As 
expected, the probability distributions of mo for the two 
subsamples, shown in Figure [TBI are consistent with each 
other. This is also explicitly demonstrated in Figure [TTl 
which shows the probability density of the difference be- 
tween the mo in the two subsamples. The difference is 
consistent with zero within la. 

We next examine subsets defined according to physical 
properties of the sources. The first criterion we apply is 
whether the source has been detected by Fermi-LAT at a 
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Fig. 16. — Probability density of mo for the subset of bright 
CGRaBS blazars: those with seconds of RA < 30 s (solid line) or 
> 30 s (dashed line). The two distributions are consistent with a 
single value. 
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Fig. 17. — Probability density of the difference between t he m ean 
modulation index mo for the two sets considered in Fig ure[H The 
difference is consistent with zero within la. 
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Fig. 18. — Probability density of mo for CGRaBS blazars in our 
monitoring sample that are (solid line) and are not (dashed line) 
included in 1LAC. The two distributions are not consistent with a 
single value. 




-0.15 



0.15 



Fig. 19. — Probability density of the difference between t he m ean 
modulation index mo for the two sets considered in Figure [T8l The 
peak of the distribution is 7a away from zero. 



significance level high enough to warrant inclusion in the 
1LAC catalog. For sources with Sq < 0.4 Jy we apply a 
cut m > m u = 0.06 and for sources with Sq > 0.4 Jy a 
cut m > mi = 0.02. The results are shown in Figures IT8l 
and[T9l The set of sources that are included in 1LAC is 
depicted by a solid line, while the set of sources that are 
not in 1LAC is depicted by a dashed line. The two are 
not consistent with each other at a confidence level of 7a 
(Figure fT9|) . with a maximum- likelihood difference of 6 
percentage points, with gamma-ray-loud blazars exhibit- 



ing, on average, a higher variability amplitude by almost 
a factor of 2 versus non gamma-ray-loud blazars. 

We also examine the variability amplitude properties 
as a function of optical spectral classification. We ana- 
lyze the subsets of CGRaBS BL Lacs and FSRQs. The 
probability densities for the mean mo of the two subsets 
are shown in Fig [201 The results for BL Lacs (FSRQs) 
are plotted as a solid (dashed) line. The two curves are 
not consistent with each other — the BL Lacs appear to 
have, on average, higher variability amplitude than the 
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Fig. 20. — Probability density of mo for BL Lac (solid line) and 
FSRQ (dashed line) CGRaBS blazars in our monitoring sample. 
The two distributions are not consistent with a single value. 
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Fig. 21. — Probability density of the difference between t he m ean 
modulation index mo for the two sets considered in Figure l20l The 
peak of the distribution (at 0.035) is more than 3<r away from zero. 



FSRQs. We verify this finding by plotting, in Figure |2"T1 
the probability density of the difference between the mo 
of BL Lacs and FSRQs. The most likely difference is 3.5 
percentage points, and it is more than 3tr away from zero. 
Note that the difference between BL Lacs and FSRQs is 
less significant than that between gamma-ray-loud and 
non gamma-ray-loud blazars. This is both because the 
most likely difference in too values between the BL Lac 
and FSRQ subsets is smaller and because the BL Lac 
sample is smaller than the gamma-ray-loud blazar sam- 
ple: only 94 BL Lacs satisfy the data cuts we impose, 




Fig. 22. — Mean m in redshift bins of 0.5 for bright (S > 0.5 Jy) 
FSRQs in our monitoring sample. 
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Fig. 23. — Probability density of mo for FSRQs in our monitoring 
sample with z < 1.0 (solid line) and z > 1.0 (dashed line) . The 
two distributions are not consistent with a single value. 



versus f90 gamma-ray-loud blazars. As a result, the 
constraints on the intrinsic distribution of modulation 
indices (i.e. on too) are stronger in the latter case. 

Finally, we examine the dependence of variability am- 
plitude on redshift. In Figure [22] we plot the mean to 
(as calculated by a simple average rather than the like- 
lihood analysis) in redshift bins of Az = 0.5 for bright 
(S > 0.4 Jy) FSRQs with known redshifts in our mon- 
itoring sample. We exclude BL Lacs from this analysis 
so as not to bias the result, as BL Lacs with known red- 
shifts are located at low z, and we have also already 
shown that they have a higher mean to compared to FS- 
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Fig. 24. — Probability density of the difference between t he m ean 
modulation index mg for the two sets considered in Figure l23l The 
peak of the distribution (at 0.023) is more than 3<r away from zero. 



RQs. Although the errors are large, there is a hint of 
a trend toward decreasing variability amplitude with in- 
creasing redshift. We further test the significance of this 
result by splitting sources in our monitored sample in 
high- and low-redshift subsets with the dividing redshift 
at z = 1 (dashed line in Figure |22|) . In the two subsets 
we also include faint (S < 0.4 Jy) sources, with the usual 
cut at m u = 0.06. The probability density for the mean 
niQ of each subset is shown in Figure [231 where the solid 
curve corresponds to low-redshift blazars and the dashed 
curve to high-redshift FSRQs. We find that low-redshift 
FSRQs have higher, on average, intrinsic modulation in- 
dices. The result is shown to be statistically significant in 
Figure [24j where we plot the probability density of the 
difference between mg in each subset. The most likely 
difference is found to be about 2.5 percentage points, 
and more than 3er away from zero. 

7. DISCUSSION AND CONCLUSIONS 

We have discussed in detail the OVRO 40-m telescope 
15 GHz monitoring program. We have presented results 
from the first two years of observations, including re- 
duced data and light curves for all sources in our moni- 
toring sample. 

We have derived the variability amplitude properties 
of all blazars in our sample through a likelihood analy- 
sis that deconvolves the intrinsic variability from scatter 
induced due to errors in individual flux density measure- 
ments and accounts for uncertainties due to finite (and 
different) sampling in each source to calculate an intrin- 
sic modulation index as well as uncertainties on its value. 
We have used these intrinsic modulation indices to study 
whether and how the variability amplitude is correlated 
with physical properties of our sources. 

We have found that the distribution of intrinsic mod- 
ulation indices is different between sources that have / 
have not been detected by Fermi in GeV gamma rays; 



between BL Lacs / FSRQs; and between FSRQs at high 
and low redshifts. 

Our most significant result is that gamma-ray-loud 
sources have a higher, on average, variability amplitude, 
as quantified by the intrinsic modulation indices, than 
non-gamma-ray-loud sources. The most likely difference 
in mean modulation index is about 6 percentage points, 
so that gamma-ray-loud sources have, on average, a vari- 
ability amplitude almost a factor of 2 higher than sources 
not found in 1LAC. The result is very significant statis- 
tically, with the maximum-likelihood difference being 7a 
away from 0. 

It is not clear whether a selection effect or an intrin- 
sic difference is responsible for this deviation between 
the two subsets. It is, for example, conceivable, that 
all CGRaBS blazars are potentially gamma-ray-loud at 
some part of their activity cycle and, given enough obser- 
vation time, all of them would enter their "flaring" state 
(that would presumably be characterized by enhanced 
broad-band luminosity, including increased flux density 
at 15 GHz) and would be detected in GeV gamma rays. If 
this is the case, then the blazars that have been detected 
by Fermi so far would be the ones that happened to have 
been in their "flaring" state during the first year of Fermi 
operations, and it would be expected that they are seen 
to have a higher, on average, variability amplitude in 
15 GHz as well. In this scenario, given more time, more 
blazars in our sample will enter at some point their "flar- 
ing" state; they will be detected in gamma rays, and the 
amplitude of their 15 GHz emission will also increase. If 
we were to repeat the same experiment after another two 
years of observation, the source numbers in the two sub- 
samples would change, but not the average population 
properties: more sources would be detected in gamma 
rays, but these sources would now also exhibit a higher 
m. The average mo of each population would not change 
appreciably, but sources would move from one category 
(non-gamma-ray-loud) to the other (gamma-ray- loud) . 
If on the other hand we have seen all blazars in our mon- 
itoring sample in all activity states, then the variability 
amplitudes of each are not expected to change appre- 
ciably if we observe them for longer, and the number 
statistics in the two categories will likely remain fixed 
(for a fixed gamma-ray flux detection threshold). In this 
scenario, the variability amplitude is the result of some 
intrinsic, persistent physical property of blazars, which 
is also related to the gamma-ray activity of the source. 

BL Lacs are found to have higher, on average, intrinsic 
modulation indices than FSRQs, by about 3.5 percent- 
age points. Due to the smaller-number statistics and 
on-average smaller difference variability, the difference 
is less significant, but still more than 3<r away from 0. 
In addition, among our FSRQ subsample, low redshift 
sources are found to have higher, on average, intrinsic 
modulation indices, with the most likely difference on 
the mean being about 2.5 percentage points, also more 
than 3cr away from 0. 

The latter difference is not easy to interpret as an indi- 
cation of source evolution, as there are competing effects 
that could affect the result in either way. On the one 
hand, sources at higher redshift have been observed for a 
shorter rest-frame time interval due to time-dilation ef- 
fects, so it is conceivable that high-redshift sources have 
not been followed through their complete activity cy- 
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cles, and their intrinsic modulation indices will increase 
as they are observed for longer. On the other hand, 
sources at higher redshift are being observed at a higher 
rest-frame frequency. Because the radio variability am- 
plitude increases with increasing frequency, this effect 
should yield higher modulation indices for higher-redshift 
sources. As our monitoring program is continued, the im- 
portance of the first effect will decrease (or, conversely, 
with long enough light-curves, we could select to look at 
shorter light curve segments for our low-redshift sources, 
corresponding to the same rest-frame time interval as for 
our highest-redshift sources) . The second effect is not af- 
fected by length of observation time, however it operates 
in the opposite direction to the observed effect. Should 
the trend persists as it is seen here (low-redshift sources 
have higher modulation indices) , this might be an indica- 
tion of source evolution with cosmic time toward higher 
variability amplitudes. 

We have found larger variability for the LAT-detected 
blazars, for the BL Lac-type blazars, and for the 
sources at lower redshift. In addition, we have noted 
the relatively modest overlap between the (EGRET- 
like) CGRaBS sources and the early LAT detections. 
T hese facts may be r elated: it has already been shown 
in lAbdo et all (|2010bf ) that the decrease in LAT effective 
area below 0.3 GeV has strongly biased the LAT detec- 
tions to the relatively hard-spectrum high-peak blazars, 
especially the BL Lacs, compared to the EGRET sam- 
ple. Indeed 1LAC contained ~ 50% BL Lacs, while 
for EGRET FSRQs outnumbered BL Lacs by > 3x. 
These BL Lacs are radio fainter and tend to be lower- 
power sources at relatively low redshift. Thus we expect 
from the higher variability amplitude found for BL Lacs 
and lower- z sources in this paper that the LAT-detected 
blazar sample should have higher average variability. We 
note however that the difference in variability amplitude 
between gamma-ray-loud and gamma-ray-quiet blazars 
in our sample is much larger than the difference between 
BL Lacs and FSRQs, so this effect cannot be attributed 
in its entirety to different BL Lac/FSRQ number ratios 
in the CGRaBS and LAT-detected blazar samples. As 
LAT exposure increases and as refinement of the event 
cuts allows more effective area at lower energy, we might 
expect an increase in high-power, high-redshift FSRQ de- 
tections, with steeper gamma-ray spectral indices. These 
sources would have variability of lower amplitude and/or 
longer observed timescale. Indeed, continued LAT ex- 
posure is detecting more CGRaBS sources and we ex- 
pect that our OVRO-monitored sample will allow an ex- 
cellent comparison of radio variability statistics on rest- 
frame timescales comparable to those now probed for the 
nearby BL Lacs. 

In conclusion, we have, for the first time, been able to 
explicitly demonstrate that the radio variability ampli- 
tude of blazars exhibits positive correlations with physi- 
cally meaningful properties of the sources. Our findings 



are important steps toward understanding the physical 
differences between blazars with otherwise similar prop- 
erties which, however, differ in their gamma-ray activity. 

The variability amplitude in radio frequencies has 
never before been considered as a differentiating property 
between blazar samples; this was largely due to practical 
purposes, as never before has such a large, preselected, 
statistically complete sample been monitored for as long 
a time and with as hi gh a cadence. The compilation of 
the CGRaBS sample (He alev et al J 120081 ) . for example, 
was based on radio flux density, radio spectral index, and 
X-ray flux; variability information was not included, not 
because it was not considered important, but rather be- 
cause such information was, at the time, unavailable. As 
a result, the CGRaBS catalog had only moderate success 
in predicting sources that would emerge as gamma-ray 
sources in the LAT era. However, by providing a pres- 
elected sample defined by robust statistical criteria, the 
CGRaBS sample has allowed us to make unprecedented 
progress in studying the population properties of blazars, 
as in this work. 

As the additional, non-CGRaBS blazars that have 
been discovered by Fermi have now been added to our 
monitored source sample, our program will allow us to 
confirm and expand these results in upcoming years. In 
addition, by establishing, thr ough the results o f this work 
as well as t h ose p resented in iPavlidou et al.l (|201lD and 
lAbdo et"aT1 ([201 If ), that there is a close connection be- 
tween gamma-ray and 15 GHz blazar emission, we are 
justified to expect that additional progress in blazar jet 
physics is to be expected through cross-correlations in the 
time domain between 15 GHz and Fermi-L AT gamma- 
ray light curves. Such cross -correlations will be di s cusse d 
in an upcoming publication M ax-Moerbeck et al.l (|2011| ). 
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